%

clc;clear; 

% Enter model parameters
    price = 1:20;                  % price of ore
    sbar  = 100;                   % initial ore stock
    delta = 0.9;                   % discount factor    

% Construct state and action spaces
    S = (0:sbar)';                 % vector of states
    X = (0:sbar)';                 % vector of actions
    n = length(S);                 % number of states
    m = length(X);                 % number of actions  

%

% Empty matrices to store values
    vs    = zeros(sbar+1, numel(price));
    xs     = zeros(sbar+1, numel(price));
    pstars = zeros(sbar+1, sbar+1, numel(price));

%

% Vectorized version
    for k=1:numel(price)
    p = price(k);
    [SS,XX] = gridmake(S,X);
    f = (p-XX./(1+SS)).*XX;
    f(XX>SS) = -inf;
    f = reshape(f,n,m);
    g = getindex(SS-XX,SS);
    g = reshape(g,n,m);
    clear SS XX
    
% Pack model data
    clear model
    model.reward     = f;
    model.transfunc  = g;
    model.discount   = delta;
    model.horizon    = inf;

% Solve model using policy iteration for each price
    [v,x,pstar] = ddpsolve(model);
    vs(:,k)=v;
    xs(:,k)=x;
    pstars(:,:,k)=pstar;
    end

%

%%%%% Plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
% Plot optimal value function
    plot(S,vs);
    title('Optimal Value Function');
    xlabel('Stock'); ylabel('Value');
    lgd1 = legend('Location',"eastoutside");
    legend("Box","off");
    title(lgd1,'data=Price','FontSize',8);
% Plot optimal policy function  
    plot(S,X(xs));
    title('Optimal Extraction Policy');
    xlabel('Stock'); ylabel('Extraction');
    lgd2=legend('Location',"eastoutside");
    legend("Box","off");
    title(lgd2,'data=Price','FontSize',8);
% Generate optimal path
    sinit = max(S); nyrs = 20;
    spaths = zeros(nyrs+1, numel(price));
    for z=1:numel(price)
    spath = ddpsimul(pstars(:,:,z),sinit,nyrs);
    spaths(:,z) = S(spath);
    end
% Plot optimal path
    plot(0:nyrs,spaths);
    title('Optimal State Path');
    xlabel('Year'); ylabel('Stock');
    lgd3=legend('Location',"eastoutside");
    legend("Box","off");
    title(lgd3,'data=Price','FontSize',8);

%
%
%
%
%
%
%
%
%